% jsinusoids
% trying to estimate power spectrum of unevenly spaced data
% to use in psd est of champ data
% latest date 19.11.2003
% uses lombscargle algorithm to estimate power spectrum
% refernce http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=993&objectType=file

function[x,y, Frq, fs] = sinu(a1,a2,a3,a4, numsamples),
A = [a1 a2 a3 a4]; %Amplitude
%A = [15 15 20 1]; %Amplitude
%fs = 48; %Sampling interval
%fs = 1/30.5;%one month fo ecco
fs = 1/5;%5 minutes of Julia Vz
%frq = fftfrq(numsamples, fs);
Frq = (fs/numsamples) * [ floor(numsamples /10 )  floor(numsamples /100 ) ];
%x = (0:179)/fs; %time axis
x = (1:numsamples)/fs;
y = A(1)*sin(2*pi*x*Frq(1)+a4*pi/4) + A(2)*sin(2*pi*x*Frq(2)+a4*pi/4) + normrnd(0,A(3),1,length(x));
%y = A(1)*sin(2*pi*x*Frq(1)+0)+ A(2)*sin(2*pi*x*Frq(2)+0);
%y = A(2)*sin(2*pi*x*Frq(2)+0);
%y = A(2)*sin(2*pi*x*Frq(2)+ pi/4) + normrnd(0,A(4),1,length(x));
% note coeff = inversion(x,y,1) just brigs back the power (test for
% inversion) 20 Oct 2005

% close
%[Txy,F] = tfestimate(ief, eef, hanning(720),1,720,1/5);

%for tiwari

%1 Generate a time series with known amplitude sinusoid with some noise
%2. use inversion to find out the amplitude and phase coef=inversion
%3. sqrt(sum(coef(2:3).^2))*sin((2*pi*x)./100 + atan(coef(3)/coef(2)));
%(will getback the shape of sinusoid best fitting to the time series.



% k = 1:1000;
% 
% p=1;
% 
% for p = 1:40:1000,
%     k(p:p+19) = 0;
% end;
% 
% p = find(k==0);
% k(p) = [];
% 
% 
% 
% % k = floor(rand([1,500])*1000);
% % 
% % p = 1;
% % for i = 1:500,
% %     for j = i+1:500,
% %         if k(i) == k(j),
% %         q(p) = i;
% %         p = p+1;
% %         disp(k(i));disp(k(j));
% %     end;
% % end;
% % end;
%  D = [t(k)' x(k)'];
%  lombscargle(D);
% 
% % 20.11.2003
% 
%load c:\manoj\download\TestData;

% Frq = [1/365 12/365];
% A = [0.5 1 1];
% 
% x = TestData(:,1);
% 
% y = A(1)*sin(2*pi*x*Frq(1)) + A(2)*sin(2*pi*x*Frq(2)) + normrnd(0,A(3),1,length(x))';
% 
% D = [x' y'];
% [px,py] = period(x,y,6,1);
% semilogx(1./px,py);
% hold on;
% plot(1./Frq,0,'r*');